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ABSTRACT 

Context. Classical T Tauri stars (CTTS) are surrounded by actively accreting disks. According to current models 
material falls along the magnetic field lines from the disk with more or less free-fall velocity onto the star, where the 
plasma heats up and generates X-rays. 

Aims. We want to quantitatively explain the observed high energy emission and measure the infall parameters from the 

data. Absolute flux measurements allow to calculate the filling factor and the mass accretion rate. 

Methods. We use a numerical model of the hot accretion spot and solve the conservation equations. 

Results. A comparison to data from XMM-Newton and Chandra shows that our model reproduces the main features 

very well. It yields for TW Hya a filling factor of 0.3 % and a mass accretion rate 2 x 10"^" Mq yr~^. 
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1. Introduction 



Methods: numerical - Stars: pre-main sequence - Stars: late-type - Stars: individual: TW Hya 



T Tauri stars are young (< 10 Myr), low mass (M* < 
3Mq), pre-main sequence stars exhibiting strong Ha emis- 
sion. The class of "Classical T Tauri stars" (CTTS) are 
surrounded by an accretion disk and are actively accret- 
ing material from the disk. The disks do not reach all the 
way to the stellar surface, rather they are truncated in the 
vicinity of the corotation radius. Infrared ( IR) observations 
typic ally yield inner radii of 0.07-0.54 AU (jMuzerolle et all 
l2003i) . consistent with the corotation radius. Disk material 
is ionised by energetic stellar radiation and - somehow - 
loaded onto the stellar magnetic fi eld lines, tr aditionally 
assumed to be dip ole-like (but see IValenti fc^ ohns-Krull 
12004 iGregorv et al.ll2006f l. Along the magnetic field accre- 
tion funnels or curtains develop a nd matter impacts onto 
the star at nea rly free-fall velocity (jUchida &: Shibatalll984l ; 
lKoenig]|[T99l[ ). This process can remove angular momen- 
tum from the star (Shu et al. 1994). Observationally the 
accr etion can be traced in the optica l in the Ha line pro - 
file (jMuzerolle et al.ll2000f ). ii i the IR liBeristain et aD l200lh 
and in the ultraviolet UV ([Herczeg et al.ll2005l ). Further 
support for this scenario comes from the measurement of 
magnetic fields in some CTTS using the technique of spec- 
tropo l arime t ry and Zeeman-broaden ing (jJohns-Krull et aTl 
119991 120031 : ISvmington et all l2005t) . The kinetic energy 
of the accretion stream is released in one or several hot 
spots close to the stellar surface, producing the observed 
veiling continuum, and also line emission in the UV and 
X-rays. The emitted UV continuum radiation was previ- 
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ously calculated bv lCalvet fc GuUbringl ([Till) and detailed 
models of the accretion geometry prove that stable states 
with tw o or more a ccretion spots on the surface can ex- 
ist (Romanova et al.l[2004 ). The UV emission has been also 
used to estimate mass accreti on rates in CTTS with typi- 
cal v alues of 1O~^M0 year~^ (jOomez de Castro fc LamzinI 
.1991; iJohns-Krull et'al.ll2000D . 

T Tauri stars are copious emitters of X-ray emission. 
Specifically, X-ray emission from q uite a few CTTS was 
detected with the Einstein ([Feigelson fc Decamplil 



119811: ■ , 

lites (Feigelson et al.l Il993t iNeuhauser et al.l 



Feigelson fc KrissI I1989D and ROSAT 



satel- 

- _ 119951 : 

iGregorio-Hetem et al.l Il998[ ). The origin of the de- 
tected X-ray emission is usually interpreted as a scaled-up 
version of coronal activity as observed for our Sun, and 
the data and their interpretation prior to the satellites 
X MM-Newton and Chandr a are s ummarised in a review 
by iFeigelson fc Montmerld (|l999t ). However, recent ob- 
servatio ns of the CTTS TW Hya (iKastner et al.l 120021 : 
Stclzcr fc SchmittI l2004f). BP Tau (ISchmitt et all |200. 



and V4046 Sgr (jCiinther et al.l l2006( l with the grating 
spectrometers onboard XMM-Newton and Chandra in- 
dicate very high plasma densities in the X-ray emitting 
regions much higher than those observed in typical coronal 
sources. This finding suggests a different origin of at least 
the soft part of the X-ray spectrum in CTTS. Simple 
estimates show that X-rays can indeed be produced i n 
the accretion spot of a typical CTTS ([Lamzi 3 I1998D . 
We present here a more detailed accretion shock model, 
which predicts individual emission lines and can thus be 
directly confronted with observations to determine model 
parameters such as the maximum plasma temperature and 
mass accretion rate. 

Unfortunately, only a few CTTS have so far been stud- 
ied in detail using high-resolution X-ray spectroscopy with 
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sufScient signal-to-noise ratio. It is therefore unknown at 
present whether the observed low forbidden to intercombi- 
nation line ratios in He-like ions as measured for the CTTS 
TW Hya, BP Tau and V4046 Sgr are typical for CTTS as a 
class. Lower resolution studies of CTTS show significant dif- 
ferences between individual stars, possibly caused by con- 
siderable individual coronal activity, and in general it is 
difhcult to disentangle cor onal and accretion contributions 
(|Robrade fc Schmitti[2006[ ) in the X-ray spectra of CTTS. 
Other suggested sources of X-rays include dense clumps 
in the stellar or disk wind of CTTS, he ated up by shock 
waves. Simulations by ' Matt et"al] (|2002t ) show sufficiently 
dense regions in the stellar magnetosphere. 

The goal of this paper is to consider the maximally pos- 
sible accretion shock contribution, to determine the physi- 
cal shock parameters, compute the emitted X-ray spectrum 
and assess to what extent other emission components are 
necessary to account for the observed X-ray spectra. The 
detailed plan of our paper is as follows: In Sect.[l]the obser- 
vations used in this study are described briefly, in Sect.[3]we 
present our model and give the main assumption and lim- 
itations, in Sect. 2] the results of the simulation are shown 
and applied to observational data, followed by a short dis- 
cussion of the main points in Sect. [5l 

2. Observations 

2.1. Stellar parameters 

With a distance of only 57 pc l|Wichmann et all 119981 ) 
TW Hya is the closest known CTTS; it is not sub- 
merged in a dark, molecular cloud like many other CTTS. 
Photometric observations show varia bility between magni- 
tude 10.9 and 11.3 in the V-band (jRucinski fc Krautteil 
[l983h . Broad Ha profiles (FWHM ^ 200 km s'^) 
were observed by iMuzerolle et al.l (|2000f ). and TW Hya 
apparently belongs to a group of objects with simi- 
lar age, the so-c alled TW Hydrae association (TWA; 
IWebb et all 119991 ). TW Hydrae's mass and radius are 



usually quoted as Mh, =0.7 Mf^, R ^, 
age as 10 Myr fron i IWebb et al" 
values are given by iBatalha et al 



1.0 
(|1999j) . 
(l2f)nl 



i?0 , and its 
Alternative 
who place 



TW Hya on the HR diagram bv iBaraffe et all (|1998[ ) 
and determine stellar parameters from the optical spec- 
trum, which fits an older (30 Myr) and smaller star 
(•-R^=0.8 i?0 ) The spectral type o f TW Hya is K7 V-Ml V 
(jWebb et al.l 119991 : IBatalha et all 120021) and the system 
is seen nearly p ole-on (Kast iier et al.l 119971 : IWilner et al.l 
I2OOOI : lAlencar fc Batalha 2G03). Moreover, TW Hya dis- 
plays variations in line profiles and veiling, which have 
been interpreted as signa t ures of accret i on sp ot rotation 
(jAlencar fc Batalhal [2003 IBatalha et alll2002f ). TW Hya 
has been obser ved in the UV with lUE, FUSE and 
HST/STIS (He rczeg et al.ll2002l ). revealing a weaUh of H2 
emission lines, consistent with the origin in the surface of a 
irradia ted disk, and in X-rays with ROSAT bylCosta et al.] 
(|2000t ). Chandra /YLKTGS ( KastnereTaDllOOl) and XMM- 
Newton/RGS (Stelz er fc SchmittI 12004'). where the grating 
data indictes a significant accretion shock contribution. 



2.2. X-ray observations 

We use high resolution spectra obtained with the Chandra 
and XMM-Newton. TW Hya was observed for 48 ks with 



the Chandra HETGS on July 18, 2000 (Chandra Obsid 
5). 'Kas tner et all (|2002f ) report atypically high densities 
measured from the Ovil and Neix triplets, and a very 
high neon abundance as observed for many active coro- 
nal sources in combination with a low iron abundance; 
the anomalously high neon a bundance of TW H ya was 
investigated in more detail bv lDrake etall (|2005l ). A dif- 
ferent approach to assess the plasma density of the emit- 
tin g material by means o f iron line ratios was performed 
bv iNess fc Schmitti (120051 ). First-order grating spectra were 
extracted applying standard CIAO 3.2 tools, positive and 
negative orders were added up. Individual emission lines in 
the HEG and M EG spectra were analyse d with the CORA 
line fitting tool (jNess fc Wichman'3l2002D . assuming modi- 
fied Lorentzian line profiles with (3 = 2.5. A flare occurring 
in the sec ond half of the observation was already mentioned 
bv lKastner et al.1 (|2002l ): for our global fitting approach we 
excluded the flaring period to avoid contamination due to 
the probable coronal origin of the flare. 

Another X-ray spectrum was taken with XMM- 
Newton on July 9, 2001 with an exposure time of 
30 ks (obs-id 0112880201) using the RGS as prime in- 
str ument. An analys i s of t his observation was presented 
in IStelzer fc Schmitti (|2004D . We newly reduced also this 
dataset with the XMM-Newton Science Analysis System 
(SAS) software, version 6.0 and applied the standard selec- 
tion criteria. The X-r ay spectral an alysis was carried out 
using XSPEC VI 1.3 (I ArnaudI 119961 ) . and CORA for fine- 
fitting. Because the line widths are dominated by instru- 
mental broadening we keep them fixed at AA ~ 0.06 A. 
The RGS spectra cover a larger wavelength range than 
HETGS spectra and include, in addition to the He-like 
triplets of Ne and O, also the N triplet. Both observations 
show the observed helium-like triplets to be incompatible 
with the low-density limit and an emission measure anal- 
ysis indicates the presence of plasma with a few million 
degrees emitting in the soft X-ray region (jKastner et al.l 
I2OO2I: IStelzer fc Schmitti I2004D . 

3. The model 

In the currently accepted accretion paradigm the material 
follows the magnetic field lines from the disk down to the 
surface of the star. Here we only model the base of the ac- 
cretion column, where the infalling material hits the stellar 
surface, is heated up in a shock and cools down radiatively. 
A sketch of the envisaged accretion scenario is shown in 
Fig. [TJ Calculations of the exc ess continuum prod uced in 
this reg i on we re carried out by ICalvet fc Gullbrin^ (1998. ) . 
iLamzinl (|1998( ) already computed the emerging soft X-ray 
emission from such shocks, but due to bin sizes of 50 A 
his results do not resolve individual lines, thus concealing 
much of the valuable diagnostic information. Nevertheless 
his models show that the hot spot produced by accretion 
can possibly produce not only the veiling, but also the soft 
X-ray emission. In our modelling we resolve all individ- 
ual emission lines, allowing us to use line ratios as sen- 
sitive tracers of the density and temperature in the emit- 
ting region, and determine the elemental abundances of the 
emitting regions. Additionally we explicitly consider non- 
equilibrium ionisation states and distinguish ion and elec- 
tron temperatures. 

We assume a one-dimensional, plane parallel geometry for 
the accretion column. This is reasonable because the filling 
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Fig. 1. Sketch iUustrating the structure on the accretion 
column. In the stellar atmosphere a thin standing shock 
front forms followed by a radiative cooling zone. 



factor /, expressing the ratio between total spot size Agpot 
and stellar surface A* — 47ri?^, is quite small as will be 
demonstrated a posteriori. The magnetic field is assumed 
to be perpendicular to the stellar surface and the material 
flows along the magnetic field lines. All turbulent fluxes are 
neglected and our model further assumes a stationary state. 
We use a two-fluid approximation attributing different tem- 
peratures to the atom/ion Tjon and to the electron compo- 
nents To of the plasma. However, both components move 
with the same bulk velocity v because they are strongly 
coupled by the microscopic electric fields. We ignore en- 
ergy transport by heat conduction which is justified again 
a posteriori (see Sect. 13. 4| and any radiative transport (see 
Sect. 



3.1. Calculation of shock front 

Once the infalling material impacts onto the stellar surface, 
a shock forms, defining the origin of the depth coordinate z 
(see Fig.[T]). The shock front it self is very thin, only of th e 
order of a few mean free paths (IZel'Dovich fc Raizeij|l967l ). 
thus contributing only very marginally to the total emis- 
sion. Therefore it is not necessary to numerically resolve its 
internal structure, rather it can be treated as a mathemati- 
cal discontinuity and the total change in the hydrodynamic 
variables over this dis continuity is given by the Rankine- 
Hugoniot conditions (IZel'Dovich fc Raizeil 119671 chap. 7, 
§ 15), which transform supersonic motion into subsonic 
motion in one single step. To simplify the numerical treat- 
ment we assume the direction of flow to be exactly paral- 
lel the magnetic field lines, so the Lorentz force does not 
influence the dynamics; we also expect the magnetic field 
to suppress heat co nduction. Following the treatment by 
ICalvet fc Gullbrind (|1998f) and iLamzinl (|1998D we assume 
the gas to expand into a vacuum behind the shock front. 
Because of larger viscous forces the strong shock forma- 
tion occurs only in the ionic component, while the electron 
component is first only adiabatically compressed and sub- 
sequently heated through electron- ion collisions. In order 
to calculate the state of the ionic plasma behind the shock 
front in terms of the pre-shock conditions only the fluxes of 
the conserved quantities mass, momentum and energy are 
required. Marking the state in front of the shock front by 



^^0 4 
5^ 
2po 



PqVo = 
PoVa = 

7)2 



PlVl 

Pi + 
2pi 



Pivf 



(1) 
(2) 

(3) 



where p denotes the total mass density of the gas and P its 
pressure. 

Requiring that the electric coupling between ions and 
electrons leads an adiabatic compression of the electron 
component, implies a temperature rise for the electronic 
plasma component of 



T — T 



(7-1) 



(4) 



with 7 = 5/3 denoting the adiabatic index. Because the 
time scale for heat transfer from ions to electrons is much 
larger than that of the ions passing through the shock front, 
ions and electrons leave the shock with vastly different tem- 
peratures. Numerical evaluation of the above equations be- 
hind the accretion shock front results in electron tempera- 
tures orders of magnitude lower than the ion temperature. 
Furthermore, the ions pass the shock front so fast that other 
degrees of freedom than kinetic are not excited and there- 
fore kinetic and ionisation temperatures of the ions sub- 
stantially differ. 

3.2. Structure of the post-shock region 

In the following section we compute how the originally dif- 
ferent kinetic temperatures of ions and electrons as well 
as the ionisation temperature equilibrate and calculate the 
emitted X-ray spectrum. 

3.2.1. Momentum balance 

In the post-shock region heat is transferred from the ionic 
to the electronic component, at the same time the gas ra- 
diates and cools down, so the energy of the gas is no longer 
conserved. However, the particle number flux j of ions (and 
atoms) 



J 



(5) 



is conserved, where n is the ion/atom number density; the 
electron number density is denoted by Uc. The total mo- 
mentum flux jp is conserved, since we ignore the momen- 
tum loss by radiation; it consists of the ion and the electron 
momentum as follows: 



ion + mcTlcV + Pc 
.,2 



jp = pmnnv + Pr 

= pmnJiv^ + rifcTion + mcUcV'' + ncicckT^ 
« pm^nv'^ + nfcTion + UcicckTc (me < mn) (6) 

with Pion and P^ denoting the thermodynamic pressure of 
the ion and the electron gas respectively, Tjon and Tc their 
temperatures; mpj denotes the mass of a hydrogen atom, mo 
the electron mass and p is the dimensionless atomic weight. 
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3.2.2. Energy balance 

Let us next consider the energy balance in the post-shock 
region. Both for the ions and the electrons the evolution of 
the energy per particle is described by an ordinary differ- 
ential equation (ODE). Starting from the thermodynamic 
relation 



TdE - PdV ^ dU 



(7) 



where S denotes the entropy and U the internal energy of 
the plasma per heavy particle we will derive this ODE for 
the ionic component. In eqn. [7] the quantity TdS — dQ 
denotes the heat flux through the boundaries of the sys- 
tem. The system "ions" looses heat by collisions with the 
colder electrons. Heat transfer is most efficient for the 
lighter ions and especially protons have a much larger col- 
lision cross sections than neutral hydrogen and they are 
by far the most abundant species in the plasma. To de- 
scri be the heat transfer betwee n ions and electrons we fol- 
low [ziPDovi^^RiiiiB (|l967l Chapter VII, §10), who give 
the heat flow LUei per unit volume per unit time as (in cgs- 
units) 

I -Lion -to ^1 ^rp • 

UJe^ = ^knnurio 1^ — (T m K) 

where fc is Boltzmann's constant and A is the Coulomb- 
logarithm 



A w 9.4-K l.SlnTc - O.Slnrie 



(8) 



The number density of hydrogen ions is the number den- 
sity of all heavy particles n multiplied by the abundance 
of hydrogen and the hydrogen ionisation fraction : 



3 , ^ 1 lion A 

= -kinXi^n, — (T m K). 



(9) 



Transforming Eq. ([7]) written per heavy particle and taking 
the time derivative results in 

dU dV dQ 
—, — I- = T— = —— = iOei ■ 

dt dt dt dt 

Using the stationarity condition ^ — ^ + ff^ = 
transforms this into an ODE with the dependent variable 
z, measured from the shock front inwards (see Fig. [T]); dif- 
ferentiation with respect to z will be indicated by '. Thus 

VU' + PvV = LUei ■ 

The internal energy U is in this case the thermal energy 
U = ^kT, the pressure P can be rewritten using the equa- 
tion of state for a perfect gas. The specific volume V is the 
inverse of the number density V = —. We thus obtain 



- kTif^ 



vnkTw 



-UJe 



(10) 



The derivation of a corresponding equation for the electron 
component is similar. Since the heat loss of the ion gas 
is a gain for the electrons this term enters with opposite 
sign, and an additional loss term Qcoi appears representing 
energy losses by collisions which excite or ionise a heavy 
particle that in turn radiates. 

It is convenient to write the electron number density as 



with Xf, denoting the number of electrons per 
heavy particle. 

V (^XekT^ + VXeUkTc ^-^ = UJei - QcolXeTl, (11) 

because cuei already includes the factor Xg^n by definition in 
eqn[9l Thus we are left with four independent variables (n, 
V, Tion and To); therefore, given Xe, the system of the four 
hydrodynamic equations ([5l [TUl and [TT|) is closed and can 
be solved by nu merical integra tion. 

According to'SpitzeH (|l965l Chapter 5.3) particle veloc- 
ities reach a Maxwellian distribution after a few mean free 
path lengths. Evaluating the conditions behind the shock 
front we find that after a few hundred meters such distri- 
butions are established separately for both the ions and the 
electrons. We therefore assume that both ions and electrons 
each have their own individual Maxwellian velocity distri- 
butions throughout our simulation. This allows us to define 
an effective kinetic temperature for electron-atom/ion colli- 
sions. Usually collisions are treated using the ion rest frame 
as reference frame, and we fold the kinetic velocities of ion 
and electron and write the resulting Maxwellian distribu- 
tion with the effective temperature Tcff = To + Tion J^"^ ■ 
This effective collision temperature is then used to calcu- 
late t he radiative loss term Qcoi with t he CHIANTI 4.2 
code (|Dere et al.lll997tlYoung et al.ll2003l ). Because the gas 
may be in a non-equilibrium ionisation state the tables pro- 
duced by the built-in CHIANTI radJoss procedure, which 
assumes kinetic and ionisation temperatures equilibrated, 
are not valid in this case, instead a spectrum with the cur- 
rent state of ionisation and effective collision temperature is 
calculated and integrated over all contributing wavelengths 
to determine the instantaneous radiative losses. For fitting 
purposes it is further necessary to perform the simulations 
with different abundances because lower metallicity signif- 
icantly lowers the cooling rate and therefore the extent of 
the post-shock cooling zone. In this approach it has to be 
assumed that Qcoi is constant over each time step. 



3.2.3. Microscopic physics 

The calculation of the ionisation state is completely decou- 
pled from the hydrodynamic equations given above. In each 
time step the density and the temperatures of both plasma 
components are fed from the hydrodynamic into the micro- 
scopic equations below. We assume changes in ionisation 
to occur only through collisions with electrons. Ions are 
ionised by electron collisions (bound-free) and recombine 
by electron capture (free-bound). So the number density 
of ionisations ri_>i-|-i per unit time from state i to i -I- 1 is 
proportional to the number density of ions in ionisation 
state i and the number density of electrons Uc] for conve- 
nience we leave out a superscript identifying the element in 
question here: 



Recombination is the reverse process: 

7*i — Pi — >i— l^o^i • 



(12) 



(13) 



The quantity Ri^j is the rate coefficient describing ion- 
isation for J = i + 1 and recombination for j = i — \. 
For element z Ri^a = Rz+i^z+2 = because 1 represents 
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the neutral atom, which cannot recombine any further, and 
z + 1 the completely ionised ion which cannot lose any more 
electrons. The cross section a for each process depends not 
only on the ion, but also on the relative velocity of ion and 
electron. On the one hand, the number of ions in state i 
decreases by ionisation to state i + 1 or recombination to 
state i — 1, on the other hand, it increases by ionisations 
from z — 1 to i and by recombination from i + 1. For an 
element with atomic number Z there is thus a set oi Z + 1 
equations 

at 



+Ri+l^i ni+i) 



(14) 



Through the electron number density rie the equations 
for all elements are coupled and together with the condi- 
tion of number conservation they provide a complete sys- 
tem of differential equations. This system can be simpli- 
fied considerably under the assumption that is constant 
during each time step At. Because hydrogen and helium, 
the main donors of electrons, are completely ionised is 
mainly given by the hydrodynamics. This assumption leads 
to one independent set of equations per element. Dividing 
by the number density of the element in question and using 
that the number density of ions in ionisation stage Z for 
element A is = n^^x'^ with abundance of element 
A and the ionisation fraction x^, finally leads to 

ne {Ri-i^i Xi-i — {Ri^i+i + Ri^i-i) Xi 



dt 



The rate coefficients R 



(15) 

are taken from 

iMazzotta et al.l (|l998l ) for dielectronic recombination, 
for the radiative recombination and the ionisation (col- 
lisional and auto-ionisation) rate we use a code from 
D.A. Verner, which is available in electronic form on 
the welfl. We calculate only the elements with Z=l-28, 
considering all ionisation states for each of them. The 
model is implemented in IDL (Interactive data language) 
and the ODEs are independently integrated using 'Isode', 
an adaptive stepsize algorithm, which is provided in the 
IDL distribution. 

3.3. Verification 

In order to check our computational procedures we con- 
sidered several special cases with known analytical solu- 
tions. This includes a pure hydrogen gas with different 
electron and ion temperatures to test the heat transfer 
treatment and the ionisation of hydrogen at constant tem- 
perature. We compare o ur calculations f o r the coUisional 
ionisation equilibrium to IMazzotta et al] (|1998D . who use 
the same data sources as we do without the corrections 
from I Verner fc l akovlevl (|l99Clt ). CI is the element where the 
largest differences occur, for all important elements all dif- 
ferences are marginal at best. In addition to these physi- 
cal tests we examined which spatial resolution is possible. 
We use an adaptive step size sampling regions with steep 
gradients sufficiently densely so that none of our physical 
variables changes by more than 5%. In order to keep the 
computation time at a reasonable level, we also enforced a 
minimum step size of 1 m. 

^ http:/ /www. pa. uky.edu/~verner/fortran. html 



3.4. Heat conduction 

Thermal conduction tends to smooth out temperature gra- 
dients. It is not included in our simulation and we use the 
models' temp erature gradie nts to estimate its importance. 
According to ISpitzeil (|l965l Chapter 5.5) the thermal heat 



flux Fr, 



iid IS 



-f"cond — KqJ — 

dz 



(16) 



Here kq = 2 x lO^^A^^ erg K''/^s^^cm^^ is the coefficient 
of thermal conductivity. Comparing the thermal heat flux 
according to this equation to the energy flux carried by 
the bulk motion, it never exceeds more than a few percent 
of the bulk motion energy transport in the main part of 
the cooling zone except for the lowest density cases. Small 
scale chaotic magnetic fields in the plasma are possible; 
they would be frozen in and expected to further suppress 
thermal conduction. 

3.5. Optical depth effects 

Our simulation assumes all lines to be optical thin. The con- 
tinuum opacity in the soft X-ray region is small, however, 
we need to check line opacities. The optical depth r(A) for 
a given line can be expressed as 



r(A) = / K{X)dl , 



(17) 



where I measures distance along the photon path and k(A) 
is the local absorption coefficient, which can be computed 
from the oscillator strength / of the line in question, the 
number density niow of ions in the lower state and the line 
profile function $(A): 



k(A, z) 



rricC 



fni^^iz)^{X,z) , 



(18) 



with e being the electron charge and c the speed of light. 
We approximate the line profile to follow a Gaussian dis- 
tribution law with the normalisation ^{X, z)d\ = 1 at 
all z, centred at Aq with the width AAf,(z): 



$(A,z) 



1 Ao(z)^ 
/7FcAAfc(z) 



exp 



A-Ao(z) ^ 
AAfc(z) ^ 



(19) 



Ao(2:) is the wavelength at line centre. Because the shocked 
gas is moving into the star, but decelerating, it is Doppler- 
shifted with the bulk velocity v{z) at depth z according 
to Ao(z) = Arest (1 + v{z)/c), whcrc Arcst is the rest wave- 
length. The broadening AAb(z) is in the case of purely ther- 
mal broadening 



AAfc(z) 



Ao(z) 2kT{z) 



(20) 



but turbulent broadening AAt(z) may additionally con- 
tribute. This cannot be calculated from the ID- 
hydrodynamics in our approach, it can only be included 
in an ad-hoc fashion. On the one hand, at the boundaries 
of the accretion shock region typical turbulent flows might 
reach velocities comparable to the bulk motion and thus 
significantly broaden the observed lines, on the other hand. 
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the magnetic fields presumed to be present should tend to 
suppress flows perpendicular to the field lines. For the cal- 
culation of the optical depth we chose AXt{z) = 10 km s~^. 
If the turbulent broadening is larger, the line profiles get 
wider and the optical depth decreases. All lines consid- 
ered in this study are excited from the ground state. Since 
in coUisionally-dominated plasmas almost all excited ion 
states decay relatively fast, we assume that all ions are in 
their ground state, niow is then the product of the ionisation 
fraction for the line producing ion, the abundance of that 
element and the total ion number density. It is important 
to consider that the line centre depends on depth because 
of the Doppler shift due to the bulk velocity. Photons emit- 
ted at line centre in deeper regions end up in the wings 
of the profile in higher layers. Since our simulation does 
not include radiative transfer all lines have to be checked 
for optical depth effects. The exact geometry and position 
of the accretion shock is still unkown, but of substantial 
importance for estimates the line optical depth. The ra- 
diation could either escape through the boundaries of the 
post-shock funnel perpendicular to the direction of fiow or 
through the shock and the pre-shock region. If the spatial 
extent of any single accretion funnel is small and it is lo- 
cated high up in the stellar atmosphere, the optical depth 
in the first scenario is small and the accretion contribution 
to the total stellar emission is large. We select this scenario 
in the further discussion. If, however, the shock is buried 
deep in the stellar atmosphere, the radiation can only es- 
cape through the shock and the thin pre-shock gas. In this 
case, denpending in the infall conditions, the optical depth 
of resonance lines can be large compared to unity. 



3.6. Limitations by ID 

Since our model is ID, all emitted photons can travel only 
up or down, they cannot leave the accretion region side- 
ways. To a first approximation half of the photons is emit- 
ted in either direction. The downward emitted photons will 
eventually be absorbed by the surrounding stellar atmo- 
sphere which will be heated by this radiation. The influ- 
ence of the surrounding atmosphere on the shock region is 
expected to be small, since the temperature and hence the 
energy flux from the surrounding atmosphere into the shock 
region is much lower. We expect the shock structure to be 
well represented in ID, but the size of the hot spots could 
be underpredicted because, depending on the geometrical 
extend of each spot, less than half of the emission escapes. 



3. 7. Boundary conditions and limitations of the model 

All our simulations start with a pre-shock temperature of 
20000 K and and the corresponding (stationary) e quilib- 
rium ionisation state (jArnaud fc Rothenflud ITOSSI ). This 
choice of temperature is motivated by studies of the pho- 
toion isation in th e accretion stream by the post-shock emis- 
sion (jCalvet fc G uUbring 19981) and analytical considera- 
tions about the electron heat flux (jZel'Dovich fc Raized 
fl96l . While we ignore heat conduction in the post-shock 
region, across the shock front the temperature gradient is 
of course large. The electrons have then mean free path 
lengths much larger than the shock front extent will there- 
fore heat up the inflowing gas. We tested the influence of 
different initial conditions and found that the post-shock 



zone depends only marginally on the chosen initial ioni- 
sation state. We terminate our simulations when the tem- 
perature drops below 12 000 K. Here the opacity begins to 
play an important role and the energy flux from the cen- 
tral core of the star is no longer negligible compared to the 
accretion flux. Here and just behind the shock front the 
accuracy is affected because the step size reaches its lower 
limit and rapid ionisation or recombination processes are 
not resolved. So the model is expected to be more accurate 
for ions which exist at high temperatures e.g. in the for- 
mation region of Ovii or Oviii, which are precisely those 
observed at X-ray wavelengths. 

3.8. Model grids 

The free fall velocity onto a star with mass M« and radius 

-R* is 




/2GM, 
«free = \/ = 617 



Typically CTTS have masses comparable to 



and radii between 



1.5 Rq and i?. 



the 

= 4 



(21) 

Sun 

Rq 



(jMuzeroUe et al.|[2003t ). because they have not yet finished 
their main sequence contractio n. Most CTTS have inner 
disk radii of 10-90 solar ra dii (iMuzero Uc ct al. 2003), for 
the specific case of TW Hyal isner et al.. (2Q06,) find an in- 
ner disk radius of « 12i?Q, thus the actual infall speed 
can almost reach the free-fall speed. Previous analyses 
indicate particle number densities of the infalling gas of 
about 10^^ cm~'^. We therefore calculated a grid of mod- 
els with infall velocities vq varying between 200 km s~^ 
and 600 km s^^ in steps of 25 km s~^, and infall densi- 
ties no varying between 10^° cm~'^ and 10^^ cm~^ with 13 
points equally spaced on a logarithmic density scale. For 
each model in the grid we then calculate the emissivity 
for selected li nes in the X - ray region using th e versi on 5.1 
of CHIANTI (iDere et al.l (I1997D . iLandi et al.l (120061) '). We 
start out with abundances from 'Greve sse fc Sauvall ()1998f) 
and iterate the model fits until convergence (see Sect. 14. 3. 5]) . 



4. Results 

4.1. Structure of the post-shocl< region 

Fig. [3 shows typical model temperature and density pro- 
files. In the (infinitesimally thin) shock, defined at depth 
km, the ion temperature suddenly rises and cools down 
directly behind the shock front because the ion gas trans- 
fers energy to the electrons. After a few kilometers both 
ions and electrons have almost identical temperatures and 
henceforth there is essentially only radiative heat loss. This 
region we refer to as the post-shock cooling zone, where 
most of the X-ray emission originates (see Fig. [T]) . During 
radiative cooling the density rises and because of momen- 
tum conservation the gas slows down at the same time 
(Eq. [5]). As more and more energy is lost from the sys- 
tem, the density and the energy loss rate increase and the 
plasma cools down very rapidly in the end. In the example 
shown in Fig. [51 the shock reaches a depth of about 950 km, 
which is much smaller than a stellar radius, thus justifying 
our simplifying assumption of a planar geometry a posteri- 
ori]. 
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Fig. 2. Temperature and density profiles for a shock calcu- 
lated with no = 10^^ cm~'^ and vq = 525 km s~^, chemical 
abundances as found in TW Hya (Sect. I4.3.3p . 



Fig. 4. fonisation state of neon for a shock calculated with 
the same starting conditions as in Fig. [51 See text for a 
detailed description. 
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Fig. 3. fon and electron temperatures in the upper post- 
shock region for a shock calculated with the same starting 
conditions as in Fig. [2l 



The region where the electron and ion temperature sub- 
stantially differ from each other is much smaller than this 
maximum depth, thus a two-fluid treatment is not strictly 
necessary for most parts of the shock. In Fig. [3] both tem- 
peratures are plotted in comparison. At the shock front the 
electrons stay relatively cool because they are only com- 
pressed adiabatically. Behind the shock front the energy 
flows from the ions to the electrons, and already at a depth 
of 5 km ions and electrons have almost identical temper- 
atures. In Fig. |4] we show the depth dependence of the 
ionisation state of neon, which produces strong lines ob- 
served by Chandra and XMM-Newton. On passing through 
the shock front neon is nearly instantaneously ionised up to 
Neix, in the following « 50 km the mean ionisation rises 
until the equilibrium is reached. The plasma then contains 
a few percent Ne xi nuclei, about 40% of the is in the form 
of Nex, the rest comes as Neix. Further away from the 
shock the ions recombine because of the general cooling of 
the plasma, following the local equilibrium closely, so the 
fraction of Nevii rises. At the maximum depth the ions 



quickly recombine to lower ionisation states. 

The maximal temperatures are proportional to Vq as can 

be seen analytically from the equation of state for a perfect 

gas: 

rp Pi "o^o 2 /oo^ 

ni ni 

These estimates are obtained using Eq. [1] [21 [21 and neglect- 
ing the initial pressure. 

A higher infall velocity leads to a deeper post-shock cool- 
ing zone since the material reaches higher temperatures 
and consequently needs longer to cool down, so it flows for 
longer times and penetrates into deeper regions. Secondly, 
lower infall densities result in shocks with a larger spatial 
extent. Since the energy losses roughly scale with the square 
of the density, a lower density will increase the cooling 
time of the gas to cool down to photospheric temperatures. 
Fig. [51 shows that cooling zone lengths between 1 km and 
10000 km can be reached depending on the chosen model 
parameters. Our model assumes a free plasma flow during 
cooling without any direct influence from the surrounding 
stellar atmosphere. As long as the inf ailing gas has a suf- 
ficiently large ram pressure, the model assumption should 
be applicable, since the surrounding atmosphere is pushed 
away and mixes with the accreted material only after cool- 
ing. For thin gases the ram pressure is lower and the gas 
needs more time for cooling down. Therefore in the deeper 
and denser layers the approximation of a freely flowing gas 
is no longer valid and the whole set of model assumptions 
breaks down. The depth, where this happens, depends on 
the stellar parameters. The shock front forms where its ram 
pressure 



(23) 



roughly equals the gas pressure of the stellar atmosphere. 
The pressure of the stellar atmosphere rises exponentially 
with depth, so, independent of the starting point, only 
shocks with small cooling length can be described by the 
hydrodynamic modelling used here. We place a cut-off at 
z = 1000 km, where the pressure of the surrounding at- 
mosphere will be larger by about an order of magnitude 
already. 
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Fig. 5. The length of the cooling zone dependent on the log- 
arithm of the infal l density lognn and t he inf all velocity uq 
(abundances from iGreyesse fc Sauvall (|1998D 1: The length 
of the cooling region is labelled. The dotted area marks 
regions where the deeper parts of the simulated shock are 
below the cut-off. 



4.2. Optical depth 

We find that in all reasonable cases the total column den- 
sity over the whole simulated region is small, in the best-fit 
solution if turns out to Nh = 10^^ cm~^, so the continuum 
opacity is small as assumed. Because emission originates at 
all depths the mean absorption column is even less. The 
line opacity along the direction of flow reaches values con- 
siderably above unity, so if the radiation passes through the 
whole post-shock region and does not escape through the 
boundaries the emssion in the resonance lines is consider- 
ably reduced (e.g. in Fig. [5]). However, as we show in the 
following, the observed line ratios can be consistently ex- 
plained in the accretion shock model, so the assumption of 
a geometry allowing most photons to escape, seems to be 
realistic. 

4.3. Application to TW Hya 

In order to model the actual X-ray data available for 
TW Hya we applied a two-step process: First we used only 
the line fluxes of selected strong lines detected in the X-ray 
spectra and determined the best-fit shock model in an at- 
tempt to assess the maximally possible shock contribution. 
Then, in a second step, we performed a global, simultane- 
ous fit to all available - lower resolution - data allowing for 
possible additional coronal contributions. 

4.3.1. Fit to line fluxes 

Ratios between emission lines of the same element allow 
a determination of the best model parameters indepen- 
dent of the elemental abundances. Specifically, the XMM- 
Newton data contain the helium-like triplets of N,0 and 
Ne, which strongly emit at plasma temperatures of a few 
MK (see Tabled]). For these three elements the correspond- 
ing Lya lines are also measured, while we do not use any 
of the Ly/3 lines because they provide relatively little addi- 
tional temperature sensitivity and are substantially weaker 



Ne IX (f + i/r) ratios 




° 200 300 400 500 600 
infall velocity Vq [km/s] 



Fig. 6. Contour are labelled with calculates Ne G-ratios. 
The observed value for TW Hya is shaded in grey (Icr- 
error) . 

than the Lya lines. For each model we therefore compute 
three line ratios for each of the elements N, O and Ne, 
i.e., the so-called R- and G-ratios defined from the He- 
like triplets (Gabriel fc Jordanlll969l : iPorquet et"alll200lh 
through R = f/i and G = [f + i)/r respectively as well 
as the ratio of the Lya to the He-like r-lines. These nine 
line ratios are compared to the data via the x^— statistics; 
the resulting contour plot of as a function of the model 
parameters ng and vq is shown in Fig. [71 Correcting for 
the absorption does n ot alter the results since N h is small 
{Nh = 3.5-1020 cm-2: lRobrade fc Schmittll2006D . The best 
model is found for the parameters vq = 525 km s~^ and 
no — 10^2 cm~^ with an unreduced — 31.9 (7 degrees 
of freedom). One has to be careful here in interpreting the 
absolute value of the because it is derived from very few 
highly significant numbers with non-Gaussian errors. 

The strong neon lines confine the fit most effectively be- 
cause their values have the smallest statistical uncertainties. 
The density is mainly restricted by R-ratios, the velocity by 
the Lya/r-ratios, which are temperature-sensitive. The Ne 
G-ratio deviates from the best fit parameters significantly 
(see Fig.ini). Its observed value is 1.1 ±0.13, which points to 
infall velocities between 300 km s^^ and 400 km s~^; fitting 
only the remaining eight ratios we obtain the same best fit 
model as before. The fit does not depend on the chosen 
background radiation temperature in the range 6000 K to 
10000 K. 

The available Chandra data include a somewhat dif- 
ferent wavelength interval. Only the triplets from O and 
Ne are covered, but the Cliandra MEG has better spectral 
resolution, allowing to reliably measure several iron lines. 
We use the three lines of Fexvii at 16.78 A, 17.05 A and 
17.09 A, and calculate the ratio of the two doublet mem- 
bers and the total doublet to the third line (see Table [T|) . 
The best fit model for this data set has vq = 575 km s~^, 
slightly higher than found for the XMM-Newton data, and 
no — 10^2 cm^'^, with an unreduced value of 8.3 for 6 
degrees of freedom. Again there is no difference between 
6000 K and 10000 K taken as temperature for a black- 
body radiation background. The neon G-ratio points again 
to lower infall velocities. Both fits have a tight lower bound 
on the density no = 10^^ cm~'^ of the infalling gas. We es- 
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Table 1. The line ratios used in the fitting process for both observations are shown with their errors. The two rightmost 
columns show results for the best-fit scenario from the simulation with no = lO^^cm"'^ for both models. 



Line ratio 


XMM/RGS 


Chandra/HETGS 


170 ~ 525 km s ^ 


vo = 575 km s ^ 


N R-ratio 


0.33 ±0.24 


n. a. 


0.00 


0.00 


N G-ratio 


0.88 ±0.31 


n. a. 


0.77 


0.75 


N Lya/N VI r 


4.04 ± 1.00 


n. a. 


2.40 


2.82 


U H-ratio 


U.UD ± U.UO 


U.U4 ± U.UO 


n no 
U.Uz 


n no 
U.Uz 


G-ratio 


0.51 ±0.07 


0.82 ± 0.22 


0.73 


0.71 


Lya/O VII r 


2.01 ±0.18 


2.19 ±0.43 


1.49 


1.97 


Ne R-ratio 


0.50 ±0.07 


0.54 ±0.08 


0.32 


0.31 


Ne G-ratio 


1.10 ±0.10 


0.94 ±0.09 


0.80 


0.75 


Ne Lya/Neix r 


0.27 ±0.04 


0.62 ±0.06 


0.26 


0.49 


Fexvii (17.09A-F-17.05A)/16.78A 


n. a. 


3.32 ±0.88 


2.25 


2.22 


Fexvii 17.09A/17.05A 


n. a. 


0.58 ±0.16 


0.81 


0.79 


red. 


4.6 


1.4 
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infall velocity Vq [km/s] 



Fig. 7. Contours are labelled with unreduced values for 
TW Hya (7 dof), data from XMM-Newton. 

timate the error as one grid point, i.e., ±25 km for vq 
and ±0.33 for log(no). The free-fall velocity for TW Hya is 
between 500 km s~^ to 550 km s^^, setting a tight upper 
bound and we adopt vq = 525 km as the best value. We 
thus conclude that the same shock model is capable of ex- 
plaining the line ratios observed in both the XMM-Newton 
RGS and Chandra HETGS grating spectra of TW Hya. 

4.3.2. Global fit 

In our second step we implement the shock emission as 
XSPEC table model and proceed with a normal XSPEC 
analysis of the medium and high resolution XMM-Newton 
spectra. The analysis of both the Chandra and XMM- 
Newton data (|Kastner et al.ll2002l : [Stelzer fc Schmittll2004[ ) 
suggested the presence of higher temperature plasma possi- 
bly from an active corona; the flare observed in the Chandra 
observation may be a signature of activity. 

To account for these additional contributions we add 
up to three thermal VAPEC models and include interstel- 
lar extinction. These additional thermal components are 
calculated in the low density limit and are meant to rep- 
resent a coronal contribution. The elemental abundances 
for all components are coupled. It would be interesting to 
check if the accreting plasma is grain depleted compared 



to the corona, but the data quality is not sufficient to 
leave more abundances as free fit parameters. Because of 
the parameter degeneracy between the emission measure of 
the cool component and the interstellar absorption column 
Nh, we kept the lat ter fixed at Nh = 3 . 5 x l O'^" cm^'^, 
a value suggested bv iRobrade fc Schmittj (|2006l where a 
detailed discussion is given) . The fit uses the data from one 
of the XMM EPIC MOS detectors, RGSl and RGS2 and 
Chandra's HEG and MEG simultaneously in the energy 
range from 0.2 keV to 10 keV. The normalisation between 
the instruments is left independent to allow for calibration 
uncertainties and possible brightness variations between the 
observations. Because of the much larger count rates the 
lower resolution MOS detectors tend to dominate the 
statistic. To balance this we include all available grating 
information, but only one (MOSl) low resolution spectrum 
in the global fit. 

The fit results for our different models are presented in 
Table [21 Model A represents the best fitting pure accretion 
shock, while models B-C include additional temperature 
components; the main improvement is obviously brought 
about by the introduction of a high temperature {kT « 
1.30 keV) component, representing the emission from a hot 
corona. "Normal" coronae usually have emission measure 
distribution s that can be d escribed by a two-temperature 
model (|Briggs fc Pvd [20031 ) and this motivates us to add 
cool low-density plasma component in model C. Although 
the reduced is only marginally smaller, we regard this 
as a better model, because it can be naturally interpreted 
in terms of a stellar corona. A third low-density component 
does clearly not improve the fit any further (model D). In 
Fig. [5] we show the recorded EPIC MOSl low resolution 
spectrum, our best-fit model and separately the accretion 
and coronal contributions. An inspection of Fig. [5] shows 
that a energies below « 1.2 keV the overall emission is 
dominated by the shock emission, while at higher energies 
the coronal contribution dominates because of the thermal 
cut-off of the shock emission. A distinction between high 
and low density plasma is possible only by examining the 
line ratios in the He-like triples or in iron lines. In broad 
band spectra cool coronal and shock plasma exhibit the 
same signatures. The ratio of their emission measures was 
therefore taken from the RGS modelling alone (model C 
and D). Our global fit reproduces the triplet ratios quite 
satisfactorily as shown in Fig. [9] for neon and in Fig. 1101 for 
the oxygen triplet. 
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Fig. 8. Data for the MOSl (black symbols) and the M0S2 
(red/grey symbols) together with our best fit model C and 
its components 
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Fig. 10. O triplet: Data from RGSl with model C (black) 
and accretion (red/grey upper line) and coronal component 
(red/grey lower line) 




13.4 13. f 

Wavelength (Angstroem) 



Fig. 9. Ne triplet: Data from RGS2 with model C (black) 
and accretion (red/grey upper line) and coronal component 
(red/grey lower line) 



The resulting infall velocity for model C is 530_4 km s ^ 
and the density tiq = (1 x 10^^ ±3 x 10^^) cm^'^ (errors are 
statistical only) . In this scenario the total flux is dominated 
by the accretion shock (3.7 x 10~^^ ergs cm~^ s~^) which 
is about four times stronger than the cool corona (1.0 x 
10^^^ erg cm~^ s""'^) and five times stronger than the hot 
corona (0.8 x lO^^^ gj.g cni-2 g-i^ 0.3-2.5 keV band. 



Table 2. Reduced values for models with zero to three 
VAPEC components with temperature kTi to kT^ in [keV] 
and one shock component. The absorption column is given 
in units of [10^° j^]. The shock model converges to vq « 



530 km s and uq 
pegs at the free-fall velocity of ; 



in all cases except A {vq 
575 km s"^). 



Model 


kTi kT2 


kTs 


Nh 


red. x'' (dof) 


A 






=3.5 


2.8 (584) 


B 




1.33 


=3.5 


1.63 (580) 


C 


0.27 


1.35 


=3.5 


1.57 (577) 


D 


0.27 0.72 


1.26 


=3.5 


1.56 (573) 



Table 3. Abu n dance of 
iGrevesse fc Sauvall (|1998l ) and 



elements relative to 
first ionisation poten- 



tials (FIP). The errors are statistical only, we estimate the 
systematic error to 15%. 



Element abundance FIP [eVJ 



c 


20+" "^ 

'-'•^"-0.03 


11.3 


N 


c;i+0 05 
"■■-"^-0.04 


14.6 


O 


25+° °^ 

'-'■■^'-'-0.01 


13.6 


Ne 




21.6 


Mg 


-0.06 


7.6 


Si 


17+0-07 

'-'■^'-0.07 


8.1 


S 


0.02" 


10.4 


Fe 


IQ+o-Oi 


7.9 



" formal 2a limit 



4.3.3. Elemental abundances 

The abundance fitting has to be performed recursively un- 
til the abundances converges, because different abundances 
lead to different cooling functions and thus change the 
whole shock structure. Specifically, we start from the set of 
eleme ntal abundances determined by iRobrade fc SchmittI 
(|2006f ) and iterate. Our global fit procedure yields abun- 
dance values (Table [3) r elative to solar abundances from 
IGrevesse fc Sauvall (|l998f ). the errors given are purely sta- 



tistical (Icr range), while we believe the systematic error to 
be about 15%. As a cross-check we compare the intensities 
of lines from different elements in our pure shock model 
from Sect. 14.3.1! and find that the abundance ratios esti- 
mated in this way roughly agree. The final abundance es- 
timates (Table [Sj show a metal depleted plasma with the 
exception of neon, which is enhanced by about a factor of 
ten compared to the other elements and nitrogen, which is 
enhanced by a factor of two. Metal de pletion has also been 
observed in the wind of TW Hya by iLamzin et al.l ([2003) 
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and was noted bv lArgiroffi et all (|2005f ) using X-ray obser- 
vations of the non-accreting quadruple system TWA 5 in 
the vicinity of TW Hya. Stelzer & Schniitt (2004) interpret 
the abundances as a sign of grain depletion, where the grain 
forming elements condensate and mainly those elements 
are accreted, which stay in the gas phase li ke the noble 
:as ne on. This is discussed in more detail by iDrake et al.l 
2005f ). who collect evidence that metal depletion can be 
also seen in the the infrared and UV, where the spectral 
distribution indicates well advanced coagulation into larger 
orbiting bodies, which may resist the inward motion of the 
accreted gas. On the other hand, stars with active corona 
often show an enhancement of elements with a high first 
ionisation potenti al (IFIP), which also le ads to an enhanced 
neon abundance (jBrinkman et al.ll200l[ ). 



4.3.4. Filling factor and mass accretion rate 

A comparison of the observed energy flux /obs (at the dis- 
tance d to the star) and the simulated flux fsim per unit area 
allows to calculate the accretion spot size Agpot through 



A 



spot 



fohsid) 

/"sim (-^* ) 



(24) 



The filling factor / is the fraction of the stellar surface 
covered by the spot: 



/ = 



A. 



/obs(rf) rf^ 
47ri?2 - /,i,,(i?,)i?2 



spot 



(25) 



and the mass accretion rate is the product of the spot size 
and the mass flux per unit area, which in turn is the product 
of the gas density po = M'tih^o and the infall velocity vq: 



dM 
~dt 



= ^spotpwo = Aspot^^^rnB_noVo 



(26) 



We assume that half of the emission is directed outward and 
can be observed, the other half is directed inwards, where 
it is absorbed. For model C the spot size is 1 x 10^" cm^, 
yielding filling factors of 0.15% and 0.3% for i?* = 1 Rq 
and i?* = 0.8 Rq respectively and accretion rates of about 
2 X 10-^° Mq yr-i. 

4.4. Shock position in the stellar atmosphere 

We now compare the ram pressure of the infalling gas to the 
stellar atmospheric pressures as calculated from PHOENIX; 
we specifically use a density profile from AMES-cond-v2.6 
with effective temperature T^s = 4000 K, surface gravity 
\ogg = 4. and solar inetalici ty (jBrott fc Hauschildtl (|2005[ ) 
based on I Allard et al.l (|200lh ). The chosen stellar parame- 
ters resemble those of typical CTTS. The shock front is 
expected to form, where the ram pressure approximately 
equals the stellar atmospheric pressure, which increases ex- 
ponentially inwards. In a strict ID-geometry photons emit- 
ted upwards out of the cooling zone will be absorbed by an 
infinite accretion column, but in a more realistic geometry 
they can pass either through the stellar atmosphere or the 
pre-shock gas as can be seen in Fig. [TJ We estimate a lower 
limit for the hydrogen column density of Nh — 10^° cm"^, 
the actually measured col umn density is 3.5 x lO'^" cm"^ 
(jRobrade fc Schmitj l2006f ) by adding the column density 
between shock front and emitting ion to the the column 



density of the pre-shock gas, which the photons penetrate 
before escaping from the accretion funnel outside the stellar 
atmosphere. The optical depth of the stellar atmosphere is 
far higher than our lower limit for the pre-shock gas. This 
estimate proves that shocks, as described by our model, are 
act ually v i sible; for a contrary view on this subject matter 
see iDrakd (|2005l) . 



5. Discussion 

The best fit parameters of our shock model to match the 
X-ray observations of TW Hya are obtained by using a 
shock with the parameters vq « 525 km s^^ and log no ~ 
12. Previously infall velocities closer to ~ 300 km s~^ 
were reported by a num ber of authors (jLamzi 3 119981 : 
ICalvet fc GullbriM |1998[ ). Other observational evidence 
also suggests lower values; in the UV emission is found in 
highly ionised emission lines (C IV, N v and O Vi) extending 
up to « 400 km s~^ and in cool ions in absorption (Fell) 
again st a hot continuum likely emerging from an accretion 
spot (jLamzin et al.l [200^ . Since the gas strongly acceler- 
ates close to the stellar surface, the density will be lower 
in the high velocity region because of particle number con- 
servation, so, depending on the geometry of the accretion 
funnel, the emission measure in this region may well be very 
small. In this case the observed lines will have weak wings 
extending to larger velocities, which are difficult to iden- 
tify observationally. We therefore regard these observations 
only as a lower bound; an upper bound is provided by the 
free-fall velocity of ^ 500 — 550 km s~^. 
Measurements of TW Hya in different wavelength re- 
gions lead to conspicuously distinct mass accretion rates. 
Generally, the published estimates far exceed the results 
of our simulation, which gives an accretion rate of « 
(2 ±0.5) X 10"^° M ( y yr- i an d filling factors of 0.2 %-0.4%. 
lAlencar fc Batalhal ((20021 ) and lBatalha et al.l (|2002f) use op- 
tical spectroscopy and photometry and state a mass accre- 
tion rate between 10~^ and 10~® Mq yr~^ and a filling fac- 
tor of a few percent. In the UV the picture is inconsistent. 
On the one hand, two e mpirical relations for lin e intensi- 
ties as accretion tracers (j Johns-Krull et all |2000[ ) indicate 
mass accretion rate s above 3 x 10~ ^ Mq yr~^ (d a ta fro m 
IValenti et ah! (|2000[ ). evaluated bv iKastner et all (|2002[ )'). 
on the other han d fitting blackbodies on the UV-veiling by 
iMuzerolle et al.l (2000. ) suggests a significantly lower value: 
4 X 10~^° M0 yr~^. A similar procedure has been ear- 
lier applied bv ICosta et al.l (|2000D with a much larg er fill- 
ing factor. The previous X-ray analyses by Kastne r et al.l 
(l2002l . M = 10-« Mq yr-i) and lStelzer fcSchmittI (pOOl 
M = 10~^^ Mq yr~^) suffer from the problem that 
they use filling factors e xtract ed from UV-me a surem ents 
(from ICosta et al.l (|2000D and IMuzerolle et~al] (|2000D re- 
spectively) and a post-shock density calculated from X-ray 
observations which does not necessarily represent the same 
region. Our simulation now is the first attempt to rely solely 
on the X-ray measurements. The different values for mass 
accretion rates are summarised in Fig. 1111 
A physical explanation for these differences goes as fol- 
lows: At longer wavelengths one observes plasma at cooler 
temperatures and in general the filling factor and mass 
accretion rates are higher, because the spot is inhomoge- 
neous with different infall velocities. Fast particles would 
be responsible for a shock region with high temperatures 
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Fig. 11. Estimated mass accretion rates sorted by the en- 
ergy of the observation. +: o ptical and UV measurern ents; 
*: X-ray, for the data from ' Stelzer fc Schmit^ (|2004[ ) the 
dotted Une connects to the UV measurement, which deliv- 
ers the filling factor, ^: this simulation; data sources: see 
text 

which is observed in X-rays, whereas in other spectral 
bands cooler areas can be detected, and therefore the 
total area and the observed total mass flux is larger. 
Accretion spots with these properties are predicted by the 
magneto-hydrodynamic simul ations recently performed by 
iRomanova et al.l (|2004l ) and iGregorv et all (|2006l ). Very 
probably some of the difference can be attributed to in- 
trinsically changing accretion rates. Simulations of the in- 
ner flow region often show highly u nstable configurations 
(jvon Rekowski fc Brandenburdl200"6l ). 

We showed that basic properties of the X-ray spectra 
from CTTS can be naturally explained by accretion on a 
hot spot, but this is not the only X-ray emission mechanism. 
To understand especially the high energy tail we had to in- 
troduce two thermal components which fit a corona as it is 
expected in late-type stars and suggested by observations of 
activity; in the case of TW Hya the shock can dominate the 
overall X-ray emission. Forthcoming high-resolution obser- 
vations will hopefully allow to extend the sample of CTTS, 
where a similar analysis is feasible. 
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